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A LOCALLY GRADIENT-PRESERVING REINITIALIZATION FOR 

LEVEL SET FUNCTIONS 


LEI LI*, XIAOQIAN XU*, AND SAVERIO E. SPAGNOLIE* 

Abstract. The level set method commonly requires a reinitialization of the level set function 
due to interface motion and deformation. We extend the traditional technique for reinitializing the 
level set function to a method that preserves the interface gradient. The gradient of the level set 
function represents the stretching of the interface, which is of critical importance in many physical 
applications. The proposed locally gradient-preserving reinitialization (LGPR) method involves the 
solution of three PDEs of Hamilton-Jacobi type in succession; first the signed distance function 
is found using a traditional reinitialization technique, then the interface gradient is extended into 
the domain by a transport equation, and finally the new level set function is achieved with the 
solution of a generalized reinitialization equation. We prove the well-posedness of the Hamilton- 
Jacobi equations, with possibly discontinuous Hamiltonians, and propose numerical schemes for 
their solutions. A subcell resolution technique is used in the numerical solution of the transport 
equation to extend data away from the interface directly with high accuracy. The reinitialization 
technique is computationally inexpensive if the PDEs are solved only in a small band surrounding 
the interface. As an important application, the LGPR method will enable the application of the 
local level set approach to the Eulerian Immersed boundary method. 

Key words. Level set function; Reinitialization; Interface gradient; Eulerian Immersed Bound¬ 
ary Method; Discontinuous Hamiltonian 


1. Introduction. The level set method [25, 24] is a classical framework used to 
accurately and elegantly evolve Lagrangian interfaces over a fixed Eulerian grid. It 
has seen very wide application in numerous fields, from fluid-structure interactions 
(e.g., lipid vesicles [30], bubbles [3], two-phase flows [36]) to image processing [22], 
computational geometry [32], computer vision [32], and materials science [19, 32]. The 
level set method involves the tracking of a level set function (j), a continuous function 
with the property that its zero level set E = {a: : (j){x) = 0} represents the Lagrangian 
interface (for instance, the boundary between two fluid phases or an immersed elastic 
structure). However, if the interface is deformed by a velocity field, for instance, then 
the gradient of the associated level set function, V((), may grow unbounded in the 
process. To reduce the associated numerical error the level set function is commonly 
reinitialized. Even if the boundary is not highly deformed, when a local level set 
method [1, 28] is applied to reduce computational costs, reinitialization is required if 
the interface encroaches the boundary of the thin computational tube. 

For many applications, only the position and curvature of the interface are 
needed, and the level set function (/) after each reinitialization may be chosen to 
be a signed distance function [36, 3, 22]. For example, in the simulation of elastic 
structures immersed in a fluid, if the tension is assumed constant (see [3]) then the 
force depends only on the curvature of the interface so that the signed distance func¬ 
tion contains sufficient information. However, in the Eulerian immersed boundary 
method [5, 6], |V(()|r represents the stretching of the elastic structure. Consequently, 
the elastic forces depend on |V0| at the interface and the signed distance function 
cannot be used to compute these forces. One solution to this problem, shown by 
Cottet & Maitre, is to avoid reinitialization altogether and to instead to renormalize 
with a particular approximation of the Dirac delta function used in interface capture 
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[6]. However, there are situations in which this strategy is inadequate. For instance, 
it would not be effective in the local level set approach to the Eulerian immersed 
boundary method. 

In this paper we develop a method for reinitializing the level set function that 
locally preserves its gradient near the Lagrangian interface. The proposed locally 
gradient-preserving reinitialization (LGPR) method involves the solution of three 
Hamilton-Jacobi equations in succession; first the signed distance function is found 
using the traditional reinitialization technique, then the cost function is obtained by 
extending the interface gradient into the domain by a transport equation, and finally 
the new level set function is achieved by the solution of a generalized reinitialization 
equation with the cost function obtained in the previous step. The steady reinitializa¬ 
tion equation is an Eikonal equation with the cost function discontinuous at the cut 
locus of the interface. We show that the ’’proper” viscosity solution (to be defined) of 
the Eikonal equation exists and is unique. We also prove that the viscosity solution 
that vanishes at the interface of the reinitialization equation converges to this proper 
viscosity solution and hence it is the level set function desired. We then propose 
numerical schemes for their fast and accurate solution. As an important application, 
the LGPR method will enable the application of the local level set approach to the 
Eulerian Immersed boundary method, which may then be comparable in cost with 
the classical immersed boundary method of Peskin [29], but with improved stability. 

LGPR consists of some equations that are very similar with some well-studied 
equations in literature. Motivated by those results, some theoretic results are new 
in this paper. For example, about the Eikonal equation with boundary conditions 
imposed on the outer boundary, some results are available for some discontinuous 
cost functions. Due to the assumptions imposed in those references, the results can’t 
be applied to our case for the proof of the uniqueness. We provide a new proof for 
the uniqueness for our special Eikonal equation. The formula for the solution of the 
general reinitialization equation, as far as we know, is new, though the generalization 
from existing results is not hard. The numerical schemes are combinations of modified 
versions of some well-known methods except that we propose a new upwind scheme 
for the transport equation for extending quantities out from the interface. 

The paper is organized as follows. In §2 we present the sequence of PDEs involved 
in locally gradient-preserving reinitialization. In §3 we show the theoretical results, 
and give explicit formulas for viscosity solutions. Numerical schemes for solving the 
equations are the topic of §4, and a few illustrative examples are provided in §6. We 
conclude with a brief summary in §7. Proofs for several claims made throughout the 
paper about the cut locus, existence and uniqueness of the proper viscosity solution of 
the Eikonal equation with a discontinuous cost function, and other issues are included 
in the appendix. 

2. Problem setup and reinitialization method. We begin by describing 
in more detail the motivation and setup of the problem, and presenting the locally 
gradient-preserving reinitialization method. For the sake of presentation, we will 
consider as a model problem a closed one-dimensional elastic interface embedded in 
, though the method could be extended into cases with several closed interfaces or 
higher dimensions without conceptual difficulty. 

Suppose 0 is a level set function such that the zero level set F agrees with the 
interface where ^ is a Lagrange coordinate and t is time. Assume that (f) > 0 

inside F and (j) < 0 outside F. In [5], it was shown that \\7(j){X{^,t),t)\/\X^{^,t)\ = 
a(^) is independent of t when (j) is convected by the velocity field, and thus if (j) 
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is constructed initially such that a = 1, |V^|r measures the tangential stretching 
(or compression) of the interface. Generically, such an elastic structure responds 
energetically to both bending and stretching deformations. The elastic force due to 
interface bending depends on the curvature k = — V • n, where n = is 

the inward-pointing normal vector at the interface, which is unchanged under any 
reinitialization scheme that preserves the location of the level set. The elastic force 
at a point x due to interface stretching, however, is given by ([5]): 

(2.1) F{x) = V (£;'(|V <^|))- V • (^i?'(|V</.|)^) V05(<^) 

= KE\\Wct>\)W(j)S{(l)) + E'\\W(l)\)n ■ WW(j> ■ (/ - nn)\Wcl)\6{(j)) 

where I is the identity operator, fin is a dyadic product and E{-) is the elastic energy 
due to stretching. The first term in (2.1) is a force due to a curved interface under a 
certain tension, while the second term is due to tension gradients along the interface. 
We introduce the stretch function 

( 2 . 2 ) xix) = \V^\ix), xGT 

defined on the interface (time dependence is ignored). Stretching occurs in regions 
where x > 1, and compression occurs where x < 1. In the above, we require two 
quantities that may be tied to the gradient of the level set function; xix) and n ■ 
VV^- {I — nn)\r = Dsx{r{s))s, where s is the arc length parameter, Dg = d/ds and 
s = r'(s) is the unit tangent vector along the surface T. 

In the process of the convection of the interface, may have become unbounded 
(usually away from the interface), or the zero level set may have drifted towards the 
boundary of a tube in the local level set method. In this situation it is necessary 
to find a new level set function that is better behaved. So as to leave the elastic 
force unchanged during this process, the stretch function x must be preserved during 
reinitialization. In theory preserving x(r(s)) is sufficient, but in numerical application 
we must also ensure that its tangential derivative is accurately preserved. We now 
formulate the reinitialization problem in a more mathematical way. 

2.1. Locally gradient-preserving reinitialization. Suppose that (/iq is a uni¬ 
formly continuous level set function, on T = {x : (j)o{x) = 0} with x G but 
not necessarily elsewhere, ipo is assumed to be positive inside the interface T and 
negative outside T. In addition we assume that T satisfies: 

Assumption I. T is a closed, nonintersecting curve which can be decomposed 
into several segments, each of which is locally analytical throughout (including at the 
segment endpoints). 

Consider an arc-length parameterization of the interface on one such segment, 
r(s) : [a,b] —>• T is locally analytical if for every sg G there is a number 

e > 0, so that the Taylor series of T about Sg converges to T in (sg — £, sg -|- e) n [a, b]. 
That the segment endpoints are also assumed to be analytical (one-sided) removes 
certain pathological behaviors [4]. The assumption on T makes physical sense for 
practical interfaces. 

We denote by U the open domain enclosed by T. The stretch function 

(2.3) X(a;) = |V<()g|(x), xGT 

is assumed to satisfy: x(r(s)) is continuous and the derivative ZIsX(r(s)) is piecewise 
continuous. We assume 0 < ci < x(r(s)) < C 2 for two constants Ci,C 2 , which is a 
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physically relevant constraint since the stretching deformation is generally bounded 
when the material is elastic. We aim to find a new level set function (j) which is 
Lipschitz continuous (the gradient is bounded), smooth (C^) in a local band around 
r, and in particular, preserves the interface gradient, |V0|(a; S F) = x(x G F). 

We are then led to the Eikonal equation, 

, . iF(a;, V(/)) = sgn(())o)(|V(^| -/(x)) = 0 a; e 

^ > c/)ix) = 0, a: e F, 


for some suitable / that has the boundary condition f{x) = x(a;) for x S F. Here 
H{x,p) = sgn((/)o(x))(|p| — f{x)) is the Hamiltonian. The sign function sgn((/)o) con¬ 
nects the level set function to the so-called viscosity solution, as will be discussed in 
the next section. While / is known on F, part of the reinitialization process will be 
first to extend / away from the interface and into the larger domain. 

In the traditional reinitialization procedure, the new level set function p is the 
signed distance function which is recovered by solving numerically a Hamilton-Jacobi 
(H-J) equation [36], 


(2.5) 


dip 

Ih 


+ sgn(^o)(|V(^| - 1) = 0, 
pix,0) = (j)oix), 


which is inadequate in our effort to preserve the interface stretch information. 

Instead, we propose continuing the process by two extra steps to find a new 
function (j) that shares its gradient with (j)o locally near F. First, we extend /(x G 
r) = x{^ G r) from the interface out into the whole domain along the characteristic 
lines of the signed distance function by a transport equation. 


( 2 . 6 ) 


— -b sgn((^)VvJ • V/ = 0, 
fix e F,r) = xix). 


The desired level set function is then obtained by solving a generalized reinitialization 
equation, 

( 2 . 7 ) +sgn(<^o)(|V(/)| -/(x)) = 0, 

(/>(x,0) = (j)oix), 

For the remainder of the paper, references to the “reinitialization equation” are to 

(2.7) ; the earlier equation, (2.5), a special case of (2.7), will be referred to as the 
traditional reinitialization equation. For convenience we have abused the notation for 
/(x,t) and the steady cost function /(x), and similarly ())(x,t) and (j). Whether we 
mean the steady solution or the pseudo-time dependent solution should be clear by 
the context. We refer to (2.5)-(2.7) as the locally gradient-preserving reinitialization 
(LGPR) method. 

The LGPR method proposed above is straight-forward and there are no immedi¬ 
ately apparent complications, but it is not obvious that the solution of (2.7) converges 
to the solution of the Fikonal equation, (2.4), or even if it exists since the cost function 
/ developed with the transport equation in (2.6) may be discontinuous. However, as 
we will show in the following section, the solutions so obtained are well-defined so that 
the reinitialization method presented here may become a basis for fast, accurate local 
level set methods. We will also numerically determine how to preserve the interface 
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gradient x so that Z?sx(r(s)) can be recovered accurately. 

Remark 1. When f{x) is known in the whole domain, the Eikonal equation, 
Eq (2.4) can be solved using, for instance, the Fast Marching Method (EMM) [31] or 
the Fast Sweeping Method (FSM) [37]. This is one approach for initializing an origi¬ 
nal level set function, and could also be used as an alternative basis for reinitialization. 


3. Theoretical results. In this section we will show that the LGPR equations 
described in the previous section are well-posed. Specifically, we will show that the 
cost function f{x) is continuous outside of a closed set consisting of arcs and vertices 
and that the Eikonal equation has a unique “proper” solution (to be clarified later) 
given the function f{x) produced using the transport equation, (2.6). A formula for 
4>{x, t) is derived, which is found to converge to the “proper” solution of the Eikonal 
equation (2.4) in finite time. Thus, (2.7) is shown to be equivalent to (2.4). 

Note that there are some results about the Eikonal equation or some Hamilton- 
Jacobi equations with discontinuous Hamiltonians, which can’t applied to our case, 
as we will see later. For the Eikonal equation, one can find some results in [27, 34, 
35, 9, 11]. In these references, one can check that their assumptions don’t apply to 
the cost function generated by (2.6). In [2], the equation Ut -f sgn{uo)H{\/u) = 0 
has been studied but (2.7) doesn’t belong to this class. In our theoretical results, the 
proof for the uniqueness of the Eikonal equation is new and the formula for ([{x, r) 
as far as we know, is derived for the first time, though it is not hard to obtain this 
formula from existing results. 

3.1. The Eikonal equation. The Eikonal equation, (2.4), is a Hamilton-Jacobi 
equation with Hamiltonian F[{x,p) = sgn(0o(2;))(bl ~ /(^))j and / is called the cost 
function. We first introduce the definition of viscosity solutions (see [7, 13]), then we 
study the continuity (and regions of discontinuity) of / from its development by (2.6), 
and finally explore the associated solutions of the Eikonal equation. 

3.1.1. Viscosity solutions of the Eikonal equation. In the general setting 
of the Eikonal equation, solutions need not exist in the classical sense. Instead, so¬ 
lutions are developed in a weaker sense; specifically, a viscosity solution is defined as 
follows. 

Definition 3.1. A viscosity sub-solution (super-solution) of H{x,'S/(j)) = 0 is 
an upper semi-continuous function (a lower semi-continuous function), if for any C°° 
function (, when ([> — C has a local maximum (minimum) at xq which is an interior 
point, then F[^{x,V({xq)) < 0 {H*{x,V(({xo)) > 0). A viscosity solution is a contin¬ 
uous function that is both a sub- and super-solution. 

In this definition, H*{x,p) = limr_j.oesssup{i7(y,q)j \\{y,q) — (x,p)\\l 2 < r} is 
the sup-envelope, and iJ* is similarly defined to be the inf-envelope. 

For example, consider ju'J = l,a; G (—1,1) with T = {—1,1} (see exercises in 
[10]) whose viscosity solution is u = 1 — ]a:l for x G [—1,1]. The viscosity solution 
of —Ju'] = —l,x G (—1,1) with r = {—1,1} is M = Ja;] — 1 for x G [—1,1]. From 
this definition, we see why sgn((/)o) appears in the Eikonal equation: the viscosity 
solution of — / = 0 can only have kinks pointing up while the viscosity solution 
of / — = 0 can only have kinks pointing down. If we write the Eikonal equation 

as = /, for r that is not convex, (p that is negative outside E may have kinks 
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pointing down and this (/) is not a viscosity solution to |V0| = /. A common error in 
the numerical literature is that the signed distance function associated with a non- 
convex curve is treated as a viscosity solution of |V(^| = 1. 

It is natural to decompose the Eikonal equation, (2.4), into interior (x G U) 
and exterior (x € \ U) problems, and to piece the two solutions together. If / is 

continuous at F, which it is as we will show, then the interior and exterior solutions 
together form a viscosity solution over the entire domain (since the equation is also 
then satisfied on the interface F). We therefore focus on the interior and exterior 
problems separately. 

The interior problem has been solved by other authors for continuous / with 
inf / > 0, and existence and uniqueness have been established [7, 15, 18]. An integral 
representation of (j) is given by: 


(3.1) 


(t){x) = inf 


ly /(7(s))c^'S 7 ( 0 ) = x,7(L) G f|. 


(see [21, 27]), where is the space of absolutely continuous self-avoiding curves, s is 
the arc-length parameter, and L is the total arc-length (which depends on 7). 

In the exterior problem, however, even if / is continuous and inf / > 0, uniqueness 
is not guaranteed. For example, both (/<i(x) = jjxj — 1] and (t> 2 (x) = 1 — jxj where 
X G K are viscosity solutions for sgn(l — jxl)(l(/)'(x)l — 1) = 0. This motivates the 
following definition: 


Definition 3.2. The “proper” viscosity solution of the Eikonal equation, (2.4), 
is defined to be the pointwise limit of (fn as n ^ 00 , where 4>n is the viscosity solution 
satisfying (2.4) in the sense of Definition 3.1 in {jxj < n} with (/'nda^j = n) = 0 for 
n G Z and n > max^gp d(0, y). 


Here d(Ai, A2) = mixeEi,yeE 2 II2: — y\\ is the distance between two sets Ei and 
£ 2 - Under this definition, for continuous f > 0, (j) is the limit of a sequence of interior 
problems and is therefore uniquely determined (see the previous discussion). Such 
a, (f) is also a viscosity solution in the general sense of Definition 3.1. A limit of the 
integral representation of as n —>■ 00 reveals the viscosity solution for all x G 
with continuous / > 0, so that a representation of (j){x) for all x is given by 


(3.2) <j){x) = sgn((()o) inf I [ f{ 3 {s))ds 7(0) = x,j{L) G f|. 

76*” L Jo J 


Unfortunately, while continuous on F, it is not guaranteed that / obtained by 
the transport equation, (2.6), is continuous in the whole domain. Before proceeding 
any further, we must therefore understand the nature of / obtained using (2.6). 

3.1.2. The cost function. Due to the method for extending / into the whole 
domain using (2.6), the behavior of the cost function / is intimately linked to the 
behavior of (p. Aujol & Aubert [2] have shown that the viscosity solution of (2.5) 
that satisfies (p(x G F, r) = 0 converges to the steady solution, which is the proper 
viscosity solution (see Definition 3.2) of sgn((()o)(|— 1) = 0. By (3.2), since the 
cost function in this case is 1, is the signed distance function: 


(3.3) 


ip(x) 


d(x, F) X G U, 
-d(x,F) x€U^\U. 
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ip is therefore 1-Lipschitz continuous, and hence differentiable almost everywhere 
(a.e.). Of particular importance is the singular set of (/?, which is most conveniently 
uncovered by studying a projection to the interface: 

Definition 3.3. Px = {y G r|d(a;,r) = d(x,y)} is the nonempty projection of x 
onto r. Let A = {x\ffPx > 2} be the set of points for which the distance is achieved 
at multiple boundary points. The part of A inside of T is called the medial axis, and 
its closure A is called the cut locus. The skeleton S is the set of centers of maximal 
circles (with order defined by inclusion) inside F. 

First note that by the assumption on F we have that the distance between 
the cut locus and the interface is always positive, d(A,r) > 0. For x ^ AUF, we have 
that 

qf' _ Pt 

(3.4) yd{x,Px)= -—, 

\x — Px\ 

since Px and d{x, Px) are both differentiable at such a point, and since P{tx + (1 — 
f)Px) = Px for 0 < t < 1. Therefore, = sgn{(j)o)'^d{x, Px) is continuously 
extended to F, and thus ip is at F. Moreover, the line x — Px is a characteristic 
line of ip due to its alignment with Vd. Therefore, / is constant along the line x — Px 
by (2.6); in addition, f{x,T) is steady when r > d(a:,F) and f{x) = x(Pa:). / is 
continuous at x since Px is. f is thus continuous for AA x d A and is given by 

/(.) = x(P.). 

Having shown / to be continuous everywhere outside of the cut locus, we are 
left now to explore x G A. It is well known that ip(x) = sgn((j)Q)d{x, Px) is not 
differentiable at any point in A. Due to the importance of the result for studying 
our method, we provided an alternative proof in the appendix for the reference. Also 
note that An U C S C A n U [20]. By the assumptions on F, the curvature k of 
F exists except at possibly a finite number of points, and even at these points the 
left and right limits of the curvature exist; thus sup|n| < oo. When U is convex, 
this provides an estimate for the distance between the cut locus and the interface: 
d{A, F) > inf{l/K} (the proof can be found in the Appendix A; the chosen convention 
is that the curvature of a circle is positive). In general, however, there is no estimate 
for d(A,F). To proceed, we must further investigate the structure of the cut locus A. 
To this end, we observe the following: 

Lemma 3.4. Let x G For any y G Px, let z(t) = ty + {1 — t)x, 0 < t < 1. 
Then, Pz(t) = {y} and z(t) ^ A. 

Proof. The only nontrivial part of this claim is that z{t) ^ A \ A. Suppose 
z{t) G A \ A for some t. Then z(t) is the center of the osculating circle of F at y, 
and the circle centered at x with radius d(x,y) contains strictly a part of F inside, 
contradicting the fact that y G Px. □ 

Remark 2. Based on this lemma and the assumptions on the boundary curve, 
we are able to get another interesting geometric property of the set A: there are no 
cusps in A with zero angles. The proof is provided in the Appendix A. 

From Lemma 3.4, we have that the points in A are the terminal points of prop¬ 
agation along the characteristic lines of ip. Since the characteristics of p meet at the 
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cut locus, / may not be well defined there, so we define 


(3.5) f{x) = inf xiy), x€ A. 

y^Px 

Here notice / is extended to both the interior and exterior of the interface, which 
means we need to discuss the points both inside and outside of F. 

Let’s first concentrate on the part inside of F. One may notice that from Theorem 
6.2 and Corollary 7.1 of [4], we have S = An U is simply connected and the union of 
finitely many points and finitely many open locally analytical curves. Moreover, for 
every point on these open locally analytical curves, it has projection of size exactly 2. 

Then, let us consider the point outside of F. Suppose U is convex, then Hn[/° = 0. 
Otherwise, A H may contain curves extending to infinity in general. Let’s call 
Un = Bn n with Bn = {|a;| < n}. Then the Domain Decomposition Lemma in [4] 
tells us the closure of the skeleton Sn of Un agrees with the skeleton of H n inside, 
for instance. Bn/ 3 - This means for the points outside of F but inside of any bounded 
domain, by the result in [4] we can also get the similar characterization as the points 
inside of F. In summarize, we have 

Lemma 3.5. The set AnU is simply connected, consisting of finitely many points 
and open curves that are locally analytical, in addition, every point on those curves 
has a projection of size 2. Meanwhile, in any bounded domain the set AoU^^ consists 
of finitely many points and open curves that are locally analytical as well. 

For a point in A with a projection of size 2, function /, dehned by (3.5), has 
limiting values from both sides of the curve (for a proof see the Appendix B). We now 
have the following theorem, which will be important in investigating existence and 
uniqueness of the viscosity solutions to the Eikonal equation and the reinitialization 
equation: 

Theorem 3.6. The function f in (3.5) is bounded by ci and C2 and is continuous 
outside the cut locus {in ]R^\A), and the cut locus is well-separated from the interface, 
d(A, F) > 0. Except at finitely many points, f has a limit when x approaches A from 
one side of A. 

3.1.3. Viscosity solutions with a discontinuous cost function. We may 

now investigate the solution of the Eikonal equation when the cost function / is 
discontinuous with properties described in Theorem 3.6. The solution is the steady 
solution of the reinitialization equation and is hence the desired level set function. 
Existence has been proven in fact for a much broader class of cost functions [9]. 
Uniqueness, however, is more challenging. Uniqueness of the Eikonal equation has 
been shown for cost functions / satisfying certain conditions not applicable to the 
present case [27, 34, 35, 9, ?], so we must develop uniqueness of the solution particular 
to the cost function / of the form in Theorem 3.6. 

To begin, since / is continuous on F, we can split (2.4) into interior and exterior 
problems, as discussed in §3.1.1. We hrst consider the interior problem, 

.3 g, |V(/)| - f{x) =0, lix&U, 

' ' ' (j){x € F) = 0. 

We have here that Hf,{x,p) = (|p| — /)* = |p| — f* and H*{x,p) = (|p| — /)* = 
\p\ ~ f*i where the sup- and inf- envelopes were defined in §3.1.1. By Theorem 3.6, 


0 < Cl </*,/*< C 2 . f* is upper semi-continuous (USC) and /* is lower semi- 
continuous (LSC). If / is continuous at x, f*{x) = /*(a;) and both /*, /» are continuous 
at X. By the way we define / in (3.5), we have / = /*. It is simple to show if (f) is 
differentiable at xq, that f*{xo) < |V(()|(xo) < f*ixo) (see [10]). This implies that if 
(j) is Lipschitz (thus differentiable a.e.), then |V())| = / outside a set of measure zero. 
By approximating / by its sup- and inf-convolutions (see [9]): 

(3.7) f{x) = esssvcpy{f{y) - \y - x\^/e}, fe{x) = essinfy{/(y) + \y- x\^/e}, 

two viscosity solutions of the Eikonal equation (2.4) are found: 


(3.8) 

(3.9) 


'7(0) = ^>7(-^-) e r|, 

(l)m{x) = mf^l^ f4-j{s))ds 7(0) =x, 7 (L) e r|. 


The proof for the existence is quite routine [9]. For the convenience of the readers, 
we provided a detailed proof in the appendix. Note that /* and /* are integrable on 
every curve 7 € Also in the appendix, it is shown that (/>m and are Lipschitz 
continuous with the Lipschitz constant C 2 . These two solutions are the maximal and 
minimal viscosity solutions of the Eikonal equation [34, 35]. The viscosity solution is 
unique if (jiM = 0m- It is clear that 0 m = 0m if and only if for every point x G U, 
there’s a sequence of curves 7 „ G with total length (where 7 n( 0 ) = x and 
ln{Ln) G r) such that 


f*iln{s))ds and lim [ (/* -/*)( 7 „(s))ds = 0 . 

Jo 

The condition implies that (j)M{x) < 0m (a^); hence the two are equal. Conversely, if 
4>m{x) = (pmix), the condition is a straightforward corollary of the definition. 

The proof of the uniqueness in the literature can’t be applied to the discontinuous 
cost function / of Theorem 3.6. However, we are able to check (3.10) directly. The 
proof is provided in the appendix B. This proof for uniqueness is new and it might be 
modified to prove the uniqueness of the Eikonal equation with a class of discontinuous 
cost functions. In the exterior problem, uniqueness of the proper viscosity solution 
may be proved by first considering the finite domain U'^ H Bn with Bn = {|a;| < n}, 
where 0„ satisfying 0n(a;) = 0 , at |a;| = n using a similar proof as in the interior prob¬ 
lem and then taking n —>■ 00 . Recalling that f = f*, 'Vfe finally have the desired result: 


(3.10) 0m (ic) = lim [ 
Jo 


Theorem 3.7. The proper viscosity solution to (2.4) with the cost function ob¬ 
tained from (2.6) is unique and is given by (3.2). It is hence C 2 -Lipschitz continuous 
and in \ A. 


3.2. The reinitialization equation. Finally, we show that the reinitialization 
equation has viscosity solutions, and the solution that is zero on T is unique and 
converges to the proper viscosity solution of the Eikonal equation (2.4). Recall the 
reinitialization equation, written more generally as 


(3.11) 


du 


-f sgn('Uo)(|VM| -g{x)) = 0, 


u(a;,0) = uo{x), 
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where ug G with UC the class of uniformly continuous functions. Here u 

could be (p or 0 and g could be 1 or /, and the Hamiltonian is written as H{x,p) = 
sgn(uo(a:))(|p| — g{x)). We assume that 0 < ci < g < C 2 and that F is the zero level 
set of uq. 

For time-dependent Hamilton-Jacobi equations, the classical solutions are not 
well-defined beyond the intersection of characteristics. For some applications, the 
multi-valued solutions are important [16]; for our purpose, we need the viscosity 
solution, whose definition is similar to the one for the Eikonal equation in (3.1), with 
the only difference being the addition of a time derivative. 

Generally, if g is not continuous, we can again approximate g by g^ and g^ and 
take the limit e —0 as we did for the Eikonal equation. Therefore, we first consider 
the case where g is continuous. Motivated by the solutions of the Eikonal equation 
and the solution provided in [ 2 ] for 5 = 1 , we construct the formula of the solution, 

{ sgn(uo)inf.yg 5 f{|Mo( 7 (T))|-k /J" g(s)ds | 7 ( 0 ) = x} r < 

sgn(uo) g{s)ds \ 7 ( 0 ) = x, 'y{L) e F} t > r^,, 

where is given by 


(3.13) 


= inf ■ 


T > d(x, F) Vr > T : inf ||mo( 7 (t))| + / 9{s)ds \ 7 ( 0 ) = x| 
7G^ L Jq J 


> inf 


J g{s)ds 


I 7(0) 



It is evident that is continuous on x and Tj, < C 2 d{x, F)/ci. This formula is closely 
related to the Lax-Hopf formula and u here has the interpretation of the value function 
with cost g (see [10]). At r^,, we must have the two expressions in (3.12) being equal, for 
otherwise, the first is always strictly larger than the second for all r > d{x, F) but this 
can’t be true at r = Lopt > d{x, F), where Lopt = liminf„_,.oo Ln and is a sequence 
such that 37 „, 7 „( 0 ) = x, 7 „(L„) e F, lim„^oo fg " gilnis))ds = inf.yg<^ g{'j{s))ds. 
In addition, we find that d(x,F) < L^pt < < C 2 d(x,F)/ci. 


Remark 3. If one were to define a simpler time, 


(3.14) 



> 


inf 

7G^ 


J g{s)ds 


I 7(0) 



then Tj, might be smaller than L^pt. When u{x, t) is given by the second formula, there 
could be no paths with lengths less than r J- e„, —>■ 0 to approximate the infimum. 
For T bigger than so-defined Tx, the first expression might be smaller than the second 
one. The dynamic programming principle in the appendix then cannot be shown. 


The two expressions given in (3.12) are continuous in both x and r and they give 
the same value at r = t^. u is then continuous in both x and r. We can also see 
that it satisfies the initial and boundary conditions of the reinitialization equation. 
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In the appendix, we verify that m is a viscosity solution. From the formula, since Tx 
is bounded by C 2 d(a;, r)/ci, we see that the solution on any compact set converges to 
the proper viscosity solution of the Eikonal equation (3.2) in finite time. 

Uniqueness of the solution may be shown under the assumption that u{x G 
r, r) = 0, which can be ensured numerically. Following [2], consider 

Ou 

(3.15) —+ (|Vu|-5)=0 xGU, 

u[x, 0) = Uq{x) X G U 
u{x G F, r) = 0, 

and 

/ 9 ?/ — 

(3.16) —-{\Wu\-g) = 0 xGE.^\U, 

u{x, 0) = uo(x) X G \U 
u(x G F, r) = 0. 

The uniqueness of the solutions for these two problems have been established if g is 
continuous and bounded below by a positive number [14]. This is enough to say that 
there is at most one viscosity solution satisfying u{x G F, r) = 0. One common mistake 
in the literature is to assume that, since sgn(uo(a: G F)) = 0, that Ut{x G F ) = 0 by 
(3.11). However, this argument is inadequate in the viscosity sense, since the set of 
viscosity solutions are unchanged by a redefinition of the sign function to a value 
sgn(O) G [—1,1]- It is sensible, however, that all viscosity solutions should have 
u(x G r,T) = 0 as the characteristics flow out of F, and fortunately we can develop 
numerical schemes to ensure u{x G r,T) = 0. 

Finally, for g equal to the cost function / obtained from (2.6), it may be discon¬ 
tinuous as previously discussed. As was done for the Eikonal equation, approximating 
/ with and fg, and taking e —>■ 0, we have the maximal and minimal solutions with 
g replaced by f* and /». And as for the Eikonal equation, these two are equal and 
the solution that is zero on the interface F is unique: 

Theorem 3.8. Assume in (3.11) that either g = 1 or g = f (the cost function 
obtained using (2.6)/ The viscosity solution that satisfies u{x G F) = 0 zs unique and 
is provided in (3.12). This solution converges to the proper viscosity solution of (2.4). 


4. Numerical schemes. We have shown in theory that the LGPR method 
yields the desired level set function. We now proceed to describe numerical schemes 
for solving the PDFs introduced in §2, with a few modifications from classical methods 
[28, 23]. We also show how subcell resolution may be used to extend the interface 
gradient away from the surface with high accuracy. First we present a numerical 
scheme for solving the transport equation which involves a second-order accurate 
upwind Essentially Non-Oscillatory (ENO) scheme with subcell resolution in space 
and Gauss-Seidel iteration in time, and then we describe a method for solving the 
reinitialization equation which involves a Godunov numerical Hamiltonian scheme in 
space and again Gauss-Seidel iteration in time. 

4.1. Numerical setup. Gonsider an Eulerian grid with uniform grid size h 
upon which the interface F is overlaid. Gridpoints {xi,yj) are defined by Xi = ih and 
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Xj = jh, with i and j taking integer values. (Unlike in the theoretical part of the 
paper above, in which x and y correspond to two points in in the remainder of 
the paper x and y are coordinates, (x,y) S M^). 

In order to approximate derivatives of possibly non-smooth functions we will rely 
on ENO finite differences (see [26]). In addition, in the solution of Hamilton-Jacobi 
equations, one-sided (upwind) derivatives are commonly used to retain causality (i.e. 
information follows the characteristics). In this paper, the following one-sided second 
order ENO finite differences will be used to approximate first derivatives. 


D 


Pij — (pi—ij 


(4.1) 


Dt4>i,j = ~ ^minmod(i:)a;a;<(>ij,U’a:x(('*+ij), 


minmodja, b} = 


0 


if ab < 0, 


sgn(a) min{|a|, |6|} else, 


where the second derivative is given by the centered difference formula 


(4.2) 


D. 


xxY'f'J — j^2 




20ij + (pi-lj) ■ 


4.2. The transport equation. Recall the definitions of the stretch function 
xix) = |V(/)|(a;) on a; G r and the inward pointing unit normal vector fi = V(/)/|V(()|. 
In solving the transport equation, we aim to accurately preserve the stretch function 
as well as its tangential derivative along the interface. 


(4.3) n-\/^cj}-{I-hh)\r = Dsx{r{s))s, 

where s is the arc-length and s(s) = r'(s) is the unit vector tangent to the surface E. 
Our strategy will be to preserve the zero level set of (p (the location of the surface E) 
and the stretch function y with at least second order accuracy, and thus IlsX(E(s))s 
is formally preserved with first order accuracy. 

In the solution of the transport equation we will use a subcell resolution (SR) 
technique to obtain the cost function / (see Eq. (2.4)). In SR, the interface is deter¬ 
mined by interpolating the obtained signed distance function (p, computing the the 
gradient there, and modifying the one-sided ENO derivatives according to the inter¬ 
face [12, 23]. To illustrate the subcell resolution technique, consider as an example 
the case < 0. Letting 

(4.4) aij = h‘^mmmod{Da:xPij, Da:xPi-i,j), 

and assuming (xr,yj) G E, xy is found by quadratic interpolation at using 

the second order derivative aij/h?. The approximations of the first derivatives are 
then given by 


sp- 

dx(po(xr,yj) « + 


n 

Dx'Popj j 
Dy4>0,ij, 
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(4.5) 


dyMxr,yj) 







where , Dy are the centered differences and 


(4.6) 



{ipij + (pi-ij) - aij/-i 


(^ipij -f sgn((^2j Dij 


where 


(4.7) ^2—ij') 4c^2j(pj_i j. 

Having now obtained the cost function on the interface, /(xr, Vj) = \/ dx(j>Q + dy(j)Q, 
the left ENO derivative of / at {i,j) is modified as 

T^—f _ fij ~ I ^ij ■ r n r ^ 

fij — ^x— 2 J ’ ’ 

while at (i — 1, j) is similarly modified. 

Next, the cost function / is extended into space by solving the transport equation, 

(4.8) |^+Vd-V/ = 0, 

where Vd = sgn((p)V(p. We will denote sgn((/jyby D^dij (where are 
acting on ip not \p\ and are modified with subcell resolution near the interface) 
and we define 


(4.9) D — maxabs'(max'(0}, min{0,£)+dy}}, 

where maxabs{a,6} = (a — 6)l{|a|>|b|} + b. For the numerical gradient we then take 

( 410 ) 77 ^- = 

\/(M,)2 + (d?j,d,,)2+e2’ V(i^.d„)2 + (i4yd,,)2+e2’ 

where e is a small parameter (chosen here to be 10“^) to avoid the case that both 
and Dy are close to zero at potentially irregular points. D^dij is so chosen to ensure 
that the information propagating to (f,j) is coming points closer to the interface, 
which follows the correct characteristic directions, and also that any oscillation in / 
on the the cut locus A (see Def. 3.3), is suppressed. Using the definitions above, a 
complete spatial discretization for the transport equation is given by 

(4.11) Lrif.,) = - D+U,) , 

where a+ = maxja, 0} and a~ = — minja, 0}. 

We now turn to the discretization of time, r, which is not real time but merely 
a parameter used to relax the system to its steady state. Different timesteps kij 
are chosen for different points to ensure stability. Since 5^^ or 51^ can be small, a 
Courant-Friedrichs-Lewy (CFL) condition for convergence requires that the time step 
be small near the interface. We use the same convention as in [23], 

(4.12) h,=CuAn{h,5f,5f}, 

where C is a constant taken small enough to ensure convergence. The largest possible 
value of C depends on the cost function, and C < 1 is sufficient when / = 1. We 
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choose C = 1/2 for the numerical examples to come, which is sufficient for the cases 
studied. Further, a Gauss-Seidel iteration is employed (values are updated using the 
newest data along the chosen sweeping directions) to allow information to propagate 
long distances in some directions with one iteration. The Gauss-Seidel iteration is 
given by 

( 4 . 13 ) fij ■<r- fij + kijLxifij) 

along the four sweeping directions {i,j) = (1 : A^, 1 : TV), {i,j) = {1 : N, N : —1 : 1), 
{i,j) = (iV : —1 : 1,1 : TV), and {i,j) = {N : —I : 1, N : —1 : 1), repeatedly. 

Remark 4 . For the application of a local level set method, it may be preferable 
to use a direct time-stepping method. One possibility is the second-order TVD Runge- 
Kutta method, 

ff+^ = //) + 

( 4 - 14 ) 2 / ~\k 


Remark 5. For the transport equation (2.6), in many references Vip is dis¬ 
cretized by centered difference. Near the interface, the signed distance function is 
so that this simple treatment is convenient and sufficient. Ftowever, once f is ex¬ 
tended into the domain where the signed distance function is not smooth, the scheme 
developed in this section is expected to be more accurate. 

4.3. The reinitialization equation. We now introduce the Godunov Hamilto¬ 
nian scheme used to solve the reinitialization equation. By a similar consideration of 
causality, we use the Godunov Hamiltonian Flij (Df Uij , Df Uij , Df Uij , 11+ Uij) , where 
Sij = sgn{uoixi,yj)), and 


(4.15) Hij{a,b,c,d) 


Sij ^•\/max{o+,6 }^-|-max{c+, d — gij'^ 
Sij ^•\/max{a“, 6+}^ -I- max{c“, d+j^ — gij'^ 


if Sij > 0, 


if siij ‘P 0, 


for spatial discretization [28, 23]. (Recall u may he p oi (f> while uq = (fo.) 

It would seem to be the case that the subcell resolution (SR) technique is not 
required in order to achieve second-order accuracy with the Godunov Hamiltonian 
scheme. Without SR, the absolute error ||Vu| — g\ does indeed scale as 0{h?) (recall 
that g can be either 1 or /, see §3.2), but the interface location is generally not 
determined at the same level of accuracy, especially when the interface gradients of u 
and uq are different. This is particularly important for (2.5) because the information 
comes from the zero set of p in (2.6). The interface gradient would be determined 
only with first order accuracy and the variation in the stretch function Dsx(J'{s))s is 
then poorly captured. Thus, subcell resolution is still required to achieve second order 
accuracy for the Godunov Hamiltonian scheme. For example, when uopjUoj-ij < 0, 
we modify Df for (xi,yj) to 


u, ■ Sf- 

Dx Uzj = + -^T[mvmod{DxxUi-i^j,DxxUij). 

0 ■ ^ 

^3 


and Df for {xi-i,yj) is similarly modified. For time discretization we again use 
Gauss-Seidel iteration using the same spatially varying timestep as in (4.12). 
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The numerical Hamiltonian ensures that information propagating to (xi, yj) comes 
from values of u closer to zero. The direction of the characteristics is preserved and 
u{x S T, r) = 0 is ensured. The scheme with this numerical Hamiltonian is monotone. 
While not identical to the case of present interest, that monotone schemes of the form 
Ut + 7?(Vw) = 0, where H is continuous, have been shown to converge to the viscosity 
solution [ 8 ] is suggestive. In practice we do observe the expected convergence. 

5. Application to Eulerian Immersed boundary method. As we have 
mentioned above, our method allows the application of local level set method to 
Eulerian immersed boundary method. We’ll basically use the local level set method 
proposed in [28]. The skeleton of the algorithm is given as following: 

0. Given the immersed interface r(^), we construct the signed distance function 
f{x). Pick three positive numbers a < j3 < ^. 

1. Let T = {x : |:/?(a;)| < 7 }. In the initialization step, we set x(a; S T) = [A^j. 
Then, use LGPR to obtain (p in T. In the reinitialization step, how LGPR is used to 
get a new level set function has been explained in detail in the above sections. 

2. Solve the fluid equations (We use Navier-Stokes equations in the computation 
example later) where the fluid body force is computed by (2.1). Evolve 4> inside T with 
c{(j))u used to convect the level set functions. We use the cutoff function introduced 
in [28] 


f 1 if |y| < /3 

ciy) = { (|y|-7)^(2|?/|+7-3/3)/(7-/3)^ if/3 < [y] < 7 ■ 

[ 0 if ll/l > 7 

4. When the interface stays far enough from the boundary of the computation 
tube or |V(()| at the interfaces stays in [1 — c, 1 + c] for a given c, go to step 3; Else 
let To be the (/? — a) neighborhood of T. We reinitialize the level set function to the 
signed distance function in Tq and move to Step 1. 

For this algorithm, a lot of practical issues have been omitted since they are 
standard or straightforward but tedious to explain. However, we would like to point 
out how to use LGPR to initialize the level set function p. The method is as following: 

In a thin tube containing the interface, for every grid point {xi,yj), supposing 
z = P{xi,yj) = A(^z) is the projection, we determine the distance by computing 
d= \ {xi,yj)—z\ and the sign of of the signed distance function ip by checking {{xi, yj) — 
z) X We record lAj(^j;)| at the point {xi,yj), which is the value of |V0o,iil- 

The distance function d{x) is obtained by the Fast Sweeping Method ([37]) and the 
sign is extended in a similar manner, ip is then determined, and we evolve (2.5) to 
improve p. The values of the cost function / at the interface are interpolated from 
\^<Po,ij\- Our method may then be applied to recover p. 

6. Numerical experiments. In this section we test the LPGR method with 
a few illustrative examples. First we show that the method achieves the expected 
accuracy in a setting in which the initial level set function and interface stretching are 
specified, and we show that the elastic force is preserved through the reinitialization 
process. Next, we simulate the behavior of a stretched membrane in a fluid using 
a local level set method, and show that reinitialization is necessary and effective 
when the membrane encroaches the boundary of the computational tube. Finally, we 
show how the method could also be used to initialize level set functions from a given 
parametric curve. 
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Fig. 1. Example 1. (a) The discontinuous cost function f constructed from the specified stretch 
X = exp(0.5p) on the interface r{9) = 1 + 0.25cos(3(9). (b) The same as (a) from above. The 

cut locus (where f is discontinuous) is more readily apparent, (c) The level set function (f after 
reinitialization, (d) Contour plots of the level set function after reinitialization. The bold line is the 
zero set corresponding to F while symbols are sample positions on the interface computed analytically. 


6.1. Example 1. We first present a typical example in which the cost function 
is discontinuous in order to show that the transport equation is successfully solved 
and that the interface gradient is preserved with the desired accuracy. Let the surface 
r be parameterized in polar coordinates by r{9) = 1 + 0.25cos(3d). We take as the 
initial level set function (j)Q{x,y) = (p{x,y) exp{0.5y) where (p is the signed distance 
function relative to L. Since the interface is non-convex, characteristics of ip from the 
interface intersect both inside and outside of F. 

Using the computational domain [—1.5,1.5]^ and grid size h = 3/256, we show 
in Figs. l(a,b) two views of the cost function / produced by solving the transport 
equation. The set where / is discontinuous (the cut locus, where characteristics 
intersect) is well captured. The cost function does not oscillate near the cut locus 
using our scheme for solving the transport equation. Figs. l(c,d) show the numerical 
solution of the reinitialized level set function </> and its contours, which converges 
rapidly even though / is discontinuous. 

To test the accuracy of the numerical method, we find all the interface points 
= (xpjj/j) or = (aiijj/p) by the interpolation formula (4.6) using the data (po 
and correspondingly points p = (xT,yj) or p = (xi,yr) using p. The gradients V/iq 
and V/) are computed using (4.5) except that and are computed using the 
level set functions themselves instead of p. We define the interface location error by 
£l = max/jaip—xrl, l^p—2/r|}, error in the interface gradient by Eq = max{|V(/o(p°) — 
\7(j){p)\} and the stretching error by £s = max{||V(/(p)| — exp(0.5?/p)|}. Fig. 2(a-c) 
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Fig. 2. Convergence results for Example 1. (a) The error in the interface location decreases 
roughly as 0{h^). (b) The error in the interface gradient decays as 0{h?). (c) The error in the 
interface stretching computed using the reinitialized level set function decays as 0{h?). 


show the decay rate of these errors with the spatial step size h. The error in the com¬ 
puted interface position decays roughly as O(h^), while \/(j) (both the direction and 
norm) is accurate to 0{h‘^) as expected. Elastic forces computed using the reinitial¬ 
ized level set function are therefore expected to carry over with first order accuracy, 
which we now probe. 


Before reinitialization 


After reinitialization 




50 


0 




-50 


Fig. 3. The x-component of the elastic force Fx in Example 1 due to the complex initial 
stretching, using grid size h = 3/128, before reinitialization (a) and after reinitialization (b). The 
force computed has different relations with the curvature for y > 0 and y < 0. The traditional 
reinitialization method would naturally lose information about interface stretching, but here we see 
that the elastic force is preserved through the reinitialization process. 


Figures 3(a-b) show the ai-component of the force associated with the prescribed 
stretching in the first example, before reinitialization using 4'q (a) and after reini¬ 
tialization using (p (b). For j/ > 0, |V(^| > 1 (the interface is in tension) while 
for 2 / < 0, |V(/)| < 1 (the interface is in compression). The force is computed on 
[—1.5,1.5]^ with h = 3/128, using a Hookean elastic energy E{x) = {Kl2){x — 1)^ 
(with K = 5). Following Peskin [29], in the description of the force (2.1) that 
might be used in an immersed boundary context we use a smoothed 6 function, 
Sh = l{|0|<3;i}(l-l-cos(7r/)/(3/i)))/(6/i). The traditional reinitialization method would 
naturally lose information about interface stretching, but here we see that the elastic 
force is preserved through the reinitialization process. 

We now test the accuracy in computing both the tangent gradient of the cost 
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Cost function gradient error 


Stretching derivative error 



h h 


Fig. 4. (a) The error in the computed cost function gradient, (b) The error in the stretching 
derivative (the tension gradient in the case of Hookean elastic energy). Both are first order accurate 
as expected. 


function and the stretching derivative (the tension gradient in the case of Hookean 
elastic energy). We define T/(p°) = V/ • (/ — nn)(p°) where h(p°) = 
and correspondingly Tx(p°) = V exp(0.5?/) • (J — hh)(p°). The cost function gradient 
error £f = max{|T/(p'^) —Ty;(p'^)|} and the stretching derivative error £ 5 / = max{|n- 
VV^- [I — hn){p) — h-VV4>^ ■ (/ —nn)(p°)|}, where n{p) —V4>/\V4>\ are then defined. 
The gradients (V/, V;/) etc) are computed using (4.5), and \7\7(j){p) is computed by 
interpolating the corresponding Hessians at the two nearby points. As shown in Fig. 4, 
the quantities are computed with roughly first order accuracy as expected. 

6.2. Example 2. As a second example, we look at the application of LGPR in 
the Eulerian immersed boundary method by simulating the dynamics of a relaxing 
membrane in a fluid, where reinitialization becomes necessary since a local level set 
method is used: the fluid is described by the Navier-Stokes equations with Reynolds 
number Re = 1, which are solved using projection method of Kim and Moin [17] on 
a two-dimensional rectangular grid with no-slip boundary conditions, while the level 
set function is only constructed and updated in a tube that contains the interface. In 
its undeformed state the membrane has an arc length parameter where 0 < ^ < 27r. 
To begin the simulation the membrane is initially stretched to an ellipse given by 
X{^) = (2 cos(^), sin(^)). The system evolves under the tension generated by the curve 
and the stretching energy is given again by linear elasticity, E{x) = {K/2){x — 1)^ 
with K = 5. 

Figs. 5(a,d) show the initial level set function and its contours, with the interface 
shown as a dark black line, on the domain [—4,4]^ with grid size h = 8/100. As 
the membrane relaxes over time, the level set function evolves to the one shown in 
Figs. 5 (b,e). The level set function is no longer adequate: the interface is close to 
the boundary of the tube where the local level set method is applied, and the level 
set function has a large gradient in the computational domain (though not yet at 
the interface). At this point we perform the LGPR algorithm in another tube that 
contains the new interface. The resulting level set is shown in Figs. 5(c,f). The 
new level set function is better behaved and the sources of possible numerical error 
associated with membrane encroachment of the boundary have vanished. The x- 
component of the force before and after reinitialization is shown in Fig. 6 . The force 
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Fig. 5. Example 2. An initially stretched elastic membrane relaxes in a Navier-Stokes fluid. 
The fluid-structure interaction is computed in a small region surrounding the membrane (a local 
level set method). (a,d) The initial level set function, and its contours, with the interface shows as a 
dark black line. (b,e) The level set function and contours after the membrane has partially relaxed; 
the membrane encroaches the boundary of the computational tube containing the interface, leading 
to the development of errors. (c,f) The level set function after reinitialization; the interface is again 
well separated from the computational boundary, and the stretching information has been preserved. 


is preserved with the expected error. 
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Fig. 6. The x-component of the force from Example 2. (a) The initial force distribution, (b) 
The force after partial relaxation but before level set reinitialization (c) The force after reinitialization 
is preserved through the LGPR process. 


6.3. Example 3: Constructing the level set function. As explained in 
Section 5, LGPR is used to initialize the level set function from a given parametric 
curve and we have done this in example 2. We now focus on the behavior of LGPR in 
initializing a level set function. We take the interface in the previous example (Recall 
that we require that the level set function satisfies |V^|(A(^)) = |Aj(^)|). 

Fig. 7 shows the numerical results on the domain [—2.5, 2.5]^ with spatial grid 
size h = 5/128. The characteristic lines intersect on the medial axis of the ellipse 
and our extension scheme gives good results at the medial axis. The interface is well 
captured by our constructed level set function (Fig. 7(c)). To test the accuracy of 
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Fig. 7. (a) The cost function constructed from the initial parametrized curve from Example 2. 
(b) The constructed level set function, (c) Contours of the level set function; the bold line is the 
zero level set, or membrane interface. 


the location and the interface gradient, we plot in Fig. 8 the location error El = 
maXp{|a;p/4 + - 1|} and the gradient error Eq = maxp{||V(;()|(p) - + x^/dl}, 

where p is a point on the interface, p = (a;^, yr) or p = (xr, Vj)- The interface location 
is retained with roughly third order accuracy while the interface gradient is retained 
with second order accuracy, as desired. 


Location error 
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(b) 


1 ~^ - 
o-lO 
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1-14 . 
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Gradient error 
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Fig. 8. (a) The interface location error for the constructed level set function in Example 3 

decays as 0{h^). (b) The interface gradient error for the constructed level set function decays as 

o(h^). 


7. Conclusion. We have extended the traditional reinitialization process of a 
level set function to a locally gradient-preserving reinitialization method, which pre¬ 
serves not only the interface location, but also information about tangential interface 
stretching. We have shown in theory that the proposed method can correctly yield a 
desired level set function. In particular, we have shown that the viscosity solutions 
of the reinitialization equation converges to the unique “proper” viscosity solution of 
the Eikonal equation, which is the desired level set function. Numerical schemes are 
proposed to solve these PDEs. The subcell resolution method for the transport equa¬ 
tion is easy to implement and can extend the interface gradient to the whole domain 
with high accuracy. Numerical examples show that our method is successful, even for 
discontinuous cost functions. In applications (Eulerian immersed boundary method, 
for example), a local level set method is desirable for computational cost reductions, 
and the LGPR method can be applied in a small region containing the interface P. 
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This reinitialization process will make possible the Eulerian formulation based on a 
local level set method for simulating physical problems where the force depends on 
stretching in addition to bending. Further details on this particular application will 
be the topic of a separate paper. 
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Appendix A. Cut locus of the interface. 

We first show that the signed distance function ip is not differentiable at any 
point in A. 

Lemma A.l. Under the assumption we stated in §5, ip is not differentiable at 
any point in A. 

Proof. Since AnF = 0, to show that ip is not differentiable we show that d{x, F) is 
not differentiable at x G ADU. (The treatment of AOU^ is similar.) By the definition 
of A, we can find yi,y 2 G h, yi A ?/2 so that d{x,yi) = d{x,y 2 ) = y}(x). Suppose 
that ip is differentiable at x. We must have |V(/?(x)| = 1. Let h = V(/ 3 (x). We have 
x + eh GU for small enough e > 0. The inequality d{x + eh, yi) > ip{x + en) > (p{x) 
then gives 


(A.l) 


d{x + eh, yi) — d{x, yi) > £ + o(e). 


and the law of cosines tells us that Z{h, x — yi) — 0. But the same argument indicates 
that Z{h, X — j/ 2 ) = 0, leading to a contradiction. □ We now show that if U is convex, 
d(A,F) = d(A,F) > inf{l//c} > 0. This follows from the following lemma about the 
local structure of the curve satisfying Assumption 1. 

Lemma A.2. For any a G A and any two points y,z G P{a), if the portion of 
F between y and z is the graph of a function over the segment yz, then 3w on the 
portion such that K{w)y,{a) > l/d{a,T) where p.(o) = 1 if a is inside F and pL{a) = —1 
otherwise. 

Proof. We take the straight line yz to be the xi-axis, and assume at point z that 
the tangent line of F has a positive slope with angle 9. The projections of y and z 
to the Xi-axis are denoted by yi and zi. As the portion of F between y and z is the 
graph of a function of Xi, then so is k: k = k{xi). By integrating k between yi and 
zi (k is integrable by the assumption on the interface), we have 

(A.2) /i(a) f K{xi)dxi = 2sm{9). 

Jyi 


By basic trigonometry, zi — yi = 2d(a, F) sin(0). Hence, we can find Xi G [yi,zi] so 


that 

(A.3) 


/r(a)K(ii) > 


2 sin( 0 ) 
2d{a, F) sin( 6 >) 


1 

d{a, F) 


□ 

Although we will not use the following “no cusp” property in this paper, this is 
still an interesting property of the geometry of cut locus. This property is especially in¬ 
teresting when we discuss uniqueness of the viscosity solution to the Eikonal equation. 


Theorem A.3. For locally analytical boundary F, there is no cusp with zero 
angle in A. 

Proof. We will sketch the proof without details. Let us take Ant/ as the example. 
Suppose that there are two edges making a zero angle at x £ Ant/. Parametrize them 
with arc lengths: 7 i(s) and 72 ( 3 ), such that 7 i( 0 ) = 72 ( 0 ) = x, 7 (( 0 ) = 72 ( 0 ) =: h 
and that Vs > 0, G (0,£), 7 i(se) A 72 (se)- By Lemma 3.4, there is an opening 
region between 71 and 72 . Choose a sequence Xfc approaching x inside the region such 
that the projection of Pxk (which is one point) approaches y. Then y G P{x). By 
Lemma 3.4, Xk-Pxk doesn’t intersect with A. As a limit line segment, x-y doesn’t 
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cross 7 i ,72 and thus is tangent to the two curves, namely fi = [y — x)/\y — x\. 71 and 
72 must be on the two sides of the line xy. Consider 71 and the portion of F on that 
side, parametrized as ro(u) with ro(0) = y. Then there exists a i5 > 0 such that the 
curvature is a smooth function on [0, d] (if necessary, redefine k( 0) = lim„_,.Q+ k(u)). 
Choose Sfc —)■ O"*" such that #P 7 i(sfc) > 2. From the fact that x — y is the tangent line 
of 7 i at X, X — 7 i(sfc) and 7 i(sfc) — y are almost parallel for sufficiently large k. Then, 
for k large enough, Wk G P"fi{sk), d{x,Wk) > d{x,y) implies that Wk must fall onto 
ro(s) for 0 < s < 5. Using 7 i(sfc) = x + Skfi + 0(s\) and P{x + Skfi) = y, we know 
Vwfc S P'yi{sk) that \wk — y\ = o(s^). By Lemma A.2 and since 6 is small, there is a 
Uk such that To(ufc) is between two points in P 7 i(sfe) andK(ufe) > l/d( 7 i(sfe), F). It is 
clear that Uk = 0{s1). The contradiction is obtained by noticing that k(0) < l/d{x, y) 
and O(s^) = |k(u/j) — k( 0)| > |l/(i( 7 i(sfe), F) — l/d{x, i/)| > Csk for some C > 0. □ 

Appendix B. Solutions of Eikonal equation with discontinuous cost 
function. 

In this section, we will prove the existence and uniqueness of the viscosity solution 
to (2.4) (a.k.a.. Theorem 3.7). 

The existence part of the proof is quite standard. However, the proof of the 
uniqueness is different from the standard argument. One may find that there is one 
essential difference in our problem: in our case, the cut locus may have bifurcation 
points. Which means, there could exist point x in the singular set of cost function, 
such that for any disc B centering at a:, B will be divided by the singular set into 
more than two parts. This makes the arguments in [27, 34, 35, 9] invalid. 

On the other hand, there is another subtle difference in our proof. Although from 
Theorem A.3, there is no cusp in the singular set of the cost function, we will not use 
this. In other words, our proof doesn’t need the classical “no cusp” assumption for 
the singular set of cost function. 

B.l. Existence. We take x S 17 as an example case. We show that (j)M is a 
viscosity solution; a similar argument applies to (pm- We have the following relation¬ 
ships: 

(B.l) Cl < essmf{f{y)\y G B{x, yj(c^ - ci)e)} < /e(x) < /,(x) 

< fix) < fix) < esssup{fiy)\y G B{x, f - ci)e)} < C 2 . 

/e increases to /* while decreases to /*, and /'^(/e) is continuous. The corre¬ 
sponding solution is given as in (3.2). One then could verify the following dynamic 
programming principle: for x gU and 0 < Si < fi(x,F), 

(B.2) p^x) = ff{s))ds + p^f{si))\-fiO) = x'^. 

The proof of this principle is omitted here. For similar arguments, one could refer to 
Chapter 10 of [10]. A direct conclusion from this principle is that 

O', 0 < p"^ < C2dix,r) 

\p‘^{x) - p^iy)\ <C 2 \x-y\ 

By the Arzela-Ascoli theorem, since U is compact, there is a subsequence p’^'‘ that 
converges uniformly to p. 

Now suppose that p — C has a local maximum at xq G U. For small d > 0, 
p — p — S\x — Xo\'^ has a strict local maximum at Xg- There is a sequence of local 
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maxima Xk for 0 ®'= —C — J|x — xqP that converges to xq by the uniform convergence 
of the function sequence. Fixing iF > 0, for fc > K, we have 

(B.4) |VC(a:fe) + 25{xk - a:o)| < f^ixk) < f’^ixk)- 

Letting fc —>• oo, since is continuous, |VC(a:o)| < Sending K ^ oo, we 

obtain |VC(a;o)| < f*{xo). For a local minimum at xq, the argument is similar. Here, 
we just use the modified function (/> — ( + 5\x — xqP and the inequalities 

(B.5) |VC(a:fe) - 25{xk - a:o)| > f'‘{xk) > feA^k) > feK^^k) 

for k > K. Thus, |VC(a;o)| > f*{xo) and ^ is a viscosity solution. 

Now, fixing any <5 > 0, we can find K > 0 so that ^{x) + <5 > (x) when k > K 

for any x G U. However, By with 7 ( 0 ) = x such that cj)^'‘ + 5 > Jq f'^i"fis))ds > 
f*{'j{s))ds. (/*( 7 (-)) is the infimum of continuous functions and thus Lebesgue 
measurable on [0,r].) Meanwhile, ^{x) — 6 < (j)‘^’‘{x) < Jq f‘^’‘i"fis))ds for k large 
enough and any 7 £ with 7 ( 0 ) = x and 7 (L) G F. Now hxing 7 and taking k -G 00 , 
the dominant convergence theorem tells us that ^{x) — S < f*{"f{s))ds. Hence: 

(B. 6 ) ^(a:)=mf^|^ /*( 7 (s))ds| 7 ( 0 ) = x, 7 (L) G r| = (/)M(a:). 

The dynamic programming principle for (j)M is still true and (j)M is also bounded 
and Lipschitz continuous with the same constants. 

B.2. Uniqueness. We again fix x G C/, and we now show that (timix) = 4>m{x). 
For any £ > 0, we can find 7 G with 7 ( 0 ) = x and "f{L) G F so that f^("f(s))ds < 
(pmix) + e. 7 has no self-intersection by the definition of 

By Lemma 3.5, except at finitely many points, all the points in A belong to some 
locally analytical curves and have a projection with size 2. Noticing that 7 is an 
injection, we pick a set E that covers the finite irregular points such that the total 
length of 7 falling into E is less than e/c 2 , where C 2 is the upper bound of /*. Then 
the remaining part A \ E has the following properties: it is the disjoint union of N 
edge portions e„, 1 < n < TV; for any x G e„, we can find a ball B{x,rx) {xx > 0) 
so that every point in H n B{x, Xx) has a projection of size 2 and A n B(x, Xx) is real 
analytic. 

Step 1. We fixst show that the cost function f has limits on both sides of the edge Cn- 

For any x G e„, B{x,Xx) is divided into two subdomains Bi and B 2 by A. Let 
Xfc G Hi n [/, Xk —>■ X. The sequence Wk = Pxk has a limit point z G Px. Further in¬ 
spection reveals that z is the only limit point of Wk since ffPx = 2 and the sequences 
in B 2 give another. This means that for every sequence in Bi converging to x, the 
projections converge to z. Hence, limj,-,.^;^^^^! f{y) = x{z). If the limit function is 
/i, by the continuity of projection on one side, fi is continuous on e„. /2 may be 
similarly defined. 

Step 2. We now decompose A into sevexal paxts so that on each paxt J(f* — f*){'y{s))ds 
can be dealt with appxopxiately. 
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Let e„ be equipped with the ID Lebesgue measure m induced by the arc length, 
and let = {a; € e„ : /i(x) = f2{x)}. Clearly, f* = /* on F„ and is closed. The 
set e„ \ = {a; e e„ : /i — /2 > 0, or /i — /2 <0} is open, thus is the disjoint union of 

countable subintervals in e„. Since the sum of lengths of these subintervals is finite, 
we can find finitely many of them, say of them, such that the measure of the 
remaining is small. For these M„ intervals, we can cover the endpoints and get Mn 
new subintervals denoted as 1 < z < Hence, we can decompose e„ into G„ with 
m{Gn) < elNc2-, Fn and Let’s consider all the M = subintervals. 

3 ( 5 i ,(5 > 0 that depend on E, < n < N, such that each subinterval li 

satisfies: Ui = {x : d{x,Ii) < ( 5 i} is divided into two domains Vi\.,Vi2-, f* € C{Vi2), 
/♦ € G(f^i), and \nix£Vi2,y&Vii f*{x) — f*{y) > 5 . We have thus decomposed A into 
the union of following sets: M = subintervals; a closed set on which /* = /*; 

and a set (union of Gn and E) with measure less than 2e/c2. 

Let Gi = {s : j(s) € lij be a subset of [ 0 , L] and G = It is clear that 

f[o L]\c(‘^* ~ /*)( 7 (®))'^® < Hence, we only need to study the integral on Gi. 

Step 3 . We now obtain a local property of the portion of 7 on Gi. 

By the local smoothness of the edges, 3 q; > 0 , ^3 < (Ii/ 2 , such that if si,S2 G 
Gi, with S2 — Si < da, then the length of U between 7(51) and 7(52) is at most 
(s2 — si) + a(s2 — si)^. If 3 s e (si, S2) so that j(s) ^ Vu, then there’s an subinterval 
[53,84] C [si,S2] such that s € [53,54] and 7((s3,S4)) n I^i = 0 . Since (I3 < Si/ 2 , 
7([s 3,S4]) can’t leave ^2- Noticing /* = / = /* on 7((s3,S4)), 

/ f*ds= / f*ds> / /♦ds + (54 - 83)5 - a(s4 - 53)^02 

•Z7[S3.S4] •Z7[S3,S4] Jj 

where J is the part of U between 7(53) and 7(54). If we pick (I3 small enough, we 
could have (54 — sajd — a(s4 — 53)^02 > 0 . We replace 7([s3, 54]) with J and get a new 
curve 7, and we see that /♦ds < f f^ds. 7 may be self-intersecting. We can modify 
it as following: If J intersects with 7(0, S3), we find the infimum of Smin on [ 0 , S3] so 
that 7(5 min) ^ d. Then we piece 7[bj '^min] and the part of J starting from 7(^min) 
together. Then, by the same method we can deal with the case when J intersects 
with 7 ((s 4,L)). Such 53,54 pairs correspond to disjoint open intervals, so we can do 
this modification at most countable many times and get another curve 71 € so 
that /*ds < f*ds. Hence, without loss of generality we can assume 7 satisfies 
this property: 

If S1,S2 G Gi,\S2 - Si| < $3, then 7([si,S2]) C Vu. 

Step 4 - With the property obtained, we perturb the curve defined on Gi so that the 
integral on Gi can be treated. 


Gi is closed and consists of countable closed intervals (they are subintervals of 
[ 0 ,L], different from the intervals on the edge portion) and a nowhere dense set G' 
(G' may have positive measure). Moreover Gi = G'i is still nowhere dense since the 
extra possible points are the endpoints of the intervals. By the assumption above, we 
may write Gi = U^^G^ for some Ni G N. [inf G^,supG^] \ Gij does not contain 
any interval of length > ^3 for any j, so that 7([inf G^, sup G^]) C Vu, and 


inf 

,v2&Gij2 


|vi - z;2| > ^3 
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Now, assume K is Gij or one of the intervals in Ci, and let si = miK,Sr = sup it' 
and £1 > 0 be fixed. We can find finitely many points Sk in K and the difference 
between two consecutive points is less than lia. We shift 7([sfc, Sfc+i]) along the normal 
of li at 7(sfc) toward Vn with distance ^4. Now, we add line segments to connect the 
endpoints. The shifted curve portions and line segments are all in l^i if <^4 is small 
enough. Denote the curve so obtained by 7'. By the uniform continuity of /* in Vn^ 
we have Jy f^^ds < f^ds + Si- However, 7' may be self-intersecting. We now modify 
it following an essentially similar procedure as before: 7' consists of 7[0, s/], 7[sr,T] 
and the shifted curves together with line segments, denoted as P. Consider the first 
shifted curve portion with the line segment Pi. Suppose k > 2 is the largest number 
such that Pk intersects with Pi. We find the first point on Pi that is on P^, then 
discard the part on Pi after this point and all curve portions P; with 1 < / < fc and 
the part on Pk that is before this point. Since the portions are finite, this process 
can be terminated, resulting in P C P, P G that connects 7(sz) and j{sr)- The 
remaining work is the same as how J was modified before. Then, we find a new 
curve 72. f*ds < f*ds. By the construction, we have removed K \ {s/,Sr} 
from Ci without adding new points. Such sets are countable, we can finish this pro¬ 
cess and obtain a curve 73, so that (^^(73) consists of countably many points and 
Xyj f*ds — fifds < e/M since £1 is arbitrary. 

Hence, we are able to construct a curve 74 with 74(0) = x and 74(T) G T (since 
it has the same endpoints as 7) such that f*ds < f^ds + 2 e < f<,ds -I- 3 £ < 
4 >Tn{x) + 4 £, which verifies the condition so that 4 >Tnix) = 4 >m(,x). 

Appendix C. The solution of the reinitialization equation. 

In this section we show that the formula given in ( 3 . 12 ) is a viscosity sub-solution 
of the level-set reinitialization equation. Similar argument shows that it is also a 
super-solution is similar. 

Consider that u—(, where ( is C°°, has a local maximum at {xq, tq) for tq > 0 . We 
show that the sub-solution condition is satisfied. If tq > Tx„ , then in a neighborhood 
of {xo,tq), we have t > Tx since as Tx is continuous on x. The solution does not 
depend on r, and we also must have CT(a:o;'7)) = 0 . The sub-solution condition is 
satisfied for xq ^ T as the condition has already been verified for the Eikonal equation. 
If xo G r, [sgn(uo)(b| - g)]*\p=vaxo),x=xo < 0 is assured. 

Consider 0 < tq < Tx^. Then uo(xo) ^ 0 since for any a: G T, we have Tx = Q < tq. 
Take uo(xo) > 0 (the result for uo(xo) < 0 is similar). Let hi = min{ro, d(a:oj T)} > 0 . 
There exists 0 < h2 < hi, such that the dynamical programming principle holds for 

h < h2- 


u{l{h), TQ-h)p / 5(7(s))(is|7(0) = Xq 

Jo 

We now show this principle. Since tq < Txg, for any £ > 0 , we can find 7 with 
7(0) = Xq such that u{xo,To) + e > uo{'j{to)) + g(-j(s))ds. By the definition, 

fh° + uo{^{to)) > u{'^{h), Tq — h) whether or not tq — h > T^^k)- is thus 

shown. 

We now show the other direction. Let B = inf..y{ (7(7(s))ds|7(0) = Xq,^{L) G 
r}, and fix an arbitrary 7 G '^,7(0) = xq. We will discuss both cases where either 
To < Txq or not. 


(C.l) 


u{xq,Tq) = inf ■ 
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Consider tq < r^g. Take /12 < hi small enough so that \x — xq\ < /12, |t —ro| < ^2 
implies t < due to the continuity of Let h < /12. We can find 7i,7i(0) = 7(ft.) 
such that u(7(h), To — h)+e > g(7i(s))(is+| uo|(7i(to — h)) since >To — h. 

Connecting 7(5) : 0 < s < h and 71, we find 72. 72(7)) = 7 i(to — h). Then, 

u{xo,To) < fj” g(j 2 (s))ds + luo(72(ro))l < fg g(7(h)) + u(j(h), tq - h) + e. 

We now assume that tq = r^g (recall that we are discussing the case where 
To < Ta;g). We must have that u{xo,to) = B. Let h < hi. If tq — h > then 

371,71(0) = 7(h), 7 i(L) G T where L < tq — h such that u{^{h),To — h) + e > 
g{”fi{s))ds. Then, connecting 7(5) : 0 < s < h and 71, we get a new curve 73 
with total length h + L < tq. We then have u{xo,to) = B < g{"/3is))ds < 

Jo + u{'y{h),To — h)+e. If Tq — h < t^(/i), we can find 7i,7i(0) = 7(h) such 

that u{"f{h),To — h) + e > ^ g(^i(s))ds + |uo|(7i(ro — h)). Connecting 7(5) : 0 < 

s < h and 71, we get 72. Since tq = r^g, B < fj°g(72(s))ds + |■Uo(72(To))| by the 
dehnition of u(xo,to). The same argument as in the previous paragraph follows. 

Combining what we have, the dynamic programming principle is true. With this 
principle the sub-solution condition is easily verified: since u{xo,to) — <^(xo,to) > 
u{'j{h),To — h) — C(7(/i),To — h), we see that 


(C.2) 


C(a;o,To) < inf • 
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Cil{h),To - h) ■ 


g(7(s))ds|7(0) = xo 


which implies that 

(C.3) Cr(a;o,To) |VC|(a;o,ro) < ^(a^o)- 
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